      program bt751
      implicit real*8 (a-h,l,o-z)
      parameter (nplay=101,ndec=12,ngrid=51,npts=ngrid*ngrid)
      parameter (np=7,mp=np+1)
      parameter (ftol=1.d-10,nrep=999)
c
      external famoeb
      character*6 name(np)
      character*24 fdate
      dimension x(np),p(mp,np),y(mp),x0(np)
      dimension is(nplay,ndec)
      dimension est(np),xx(np)
      dimension bllf(nrep),boot(np,nrep),
     &          btmp(nrep),bs(np),bss(np)
      common /mats/is
      common /est/est
c
      data name/'SB2','SZ2','B1','PZ1','SB1','SZ1','A2'/
c
      open(unit=11,file='ra2011.data')   !with card sleeves on boards
      open(unit=10,file='bt751.dat')
c
      read(11,*)((is(i,n),n=1,ndec),i=1,nplay)
c
c-----------------------------------------------------------------
      write(*,'('' ***** bt751d******'',20x,a24,
     &//"    Risk & Ambiquity Experiment")')fdate()
      write(*,'(/"  All 2011 Sessions"//)')
c
      call fgrid
         call fch
c
c Set parameters
         x0(1) = dlog(1.d0/0.01d0 - 1.d0)
         x0(2) = dlog(1.d1/9.d0 - 1.d0)
         x0(3) = dlog(1.d0/0.99d0 - 1.d0)
         x0(4) = dlog(0.5d0/0.11d0 - 1.d0) 
         x0(5) = dlog(1.d0/0.26d0 - 1.d0)
         x0(6) = dlog(1.d0/0.11d0 - 1.d0)
         x0(7) = dlog(1.d0/0.29d0 - 1.d0)
c
         delta=0.1d0
         iter=2000  
         call init(x0,p,y,famoeb,delta)
         call amoeba(p,y,mp,np,np,1.d-2*ftol,famoeb,iter)
         call avg(p,x,np,mp)
         z0 = -famoeb(x)
         call m2h(x,est,np)
         write(*,'("est = ",4g13.5)')(est(k),k=1,np)
         write(*,'(/"Initial LL = ",g14.6)')z0
         write(10,910)0,0.0,z0,(est(k),k=1,np)
c 
         call savephi
c         stop
c
c begin interations
c
      dseed = 90001.d0
      do 500 it=1,nrep
   25 xseed = dseed
         call data(it,dseed)
         call fch
c
c
c ramix751 model

         delta=0.25d0
         iter=5000
         call init(x0,p,y,famoeb,delta)
         call amoeba(p,y,mp,np,np,ftol,famoeb,iter)
         if (iter.ge.5000)  then
           write(*,'(5x,"irep =",2x,i5)')it
           do 21 i=1,mp
             write(*,'(10g12.4)')(p(i,j),j=1,np),y(i)
   21      continue
           go to 25
         endif
         call avg(p,x,np,mp)
         bllf(it) = -famoeb(x)
         call m2h(x,xx,np)
         do 31 k=1,np
            boot(k,it) = xx(k)
   31    continue
         write(10,910)it,xseed,bllf(it),(boot(k,it),k=1,np)
  500 continue
  910   format(i4,2x,f20.1,10(2x,g14.6))
c
c --- process the bootstrap results ---
c
         t1=0.05d0*float(nrep)
         nlo=t1
         nlo=1+nlo/2
         nhi=nrep-nlo+1
c
      do 160 j=1,np
         bs(j)=0.d0
         bss(j)=0.d0
  160 continue
         blm=0.0d0
         blms=0.0d0
      do 162 k=1,nrep
         blm=blm+bllf(k)
         blms=blms+bllf(k)*bllf(k)
         do 162 j=1,np
            bs(j)=bs(j)+boot(j,k)
            bss(j)=bss(j)+boot(j,k)*boot(j,k)
  162 continue
      blm=blm/float(nrep)
      blms=(blms-float(nrep)*blm*blm)/float(nrep-1)
      blms=dsqrt(blms)
      do 165 j=1,np
         bs(j)=bs(j)/float(nrep)
         bss(j)=(bss(j)-float(nrep)*bs(j)*bs(j))/float(nrep-1)
         bss(j)=dsqrt(bss(j))
  165 continue
      write(*,810)nrep
      do 167 j=1,np
         tratio=est(j)/bss(j)
         do 166 k=1,nrep
  166       btmp(k)=boot(j,k)
         call sort(nrep,btmp)
         xlo=btmp(nlo)
         xhi=btmp(nhi)
         write(*,811)name(j),bs(j),bss(j),tratio,xlo,xhi
         write(*,813)est(j)
  167 continue
      write(*,'(//"The mean and std dev of log-likelihood:")')
      write(*,812)blm,blms
  810 format(//'*** BOOTSTRAP RESULTS ***  nrep=',i5//
     &  'name ','  estimate  ','  std dev   ','  t-ratio   ',
     &      2x,'---95 percent conf. int.--- ')
  811 format(a6,2x,2g13.6,2x,2g11.4,x,2g13.6)
  812 format(/2g14.6)
  813 format(7x,g14.6)
c
      end
c***********************************************************************
      subroutine m2h(x,xx,np)
      implicit real*8 (a-h,o-z)
      dimension x(np),xx(np)
c
      s0 = 0.05d0
      xx(1) = 1.d0/(1.d0+dexp(x(1)))+s0
      xx(2) = 9.95d0/(1.d0+dexp(x(2)))+s0
      xx(3) = 1.d0/(1.d0+dexp(x(3)))
      xx(4) = 0.5d0/(1.d0+dexp(x(4)))+0.5d0
      xx(5) = 1.d0/(1.d0+dexp(x(5)))+s0
      xx(6) = 1.d0/(1.d0+dexp(x(6)))+s0
      xx(7) = 1.d0/(1.d0+dexp(x(7)))
c
      return
      end
c***********************************************************************
      subroutine fgrid
      implicit real*8 (a-h,o-z)
      parameter (ndec=12,ngrid=51)
      dimension w0(ngrid),w(ngrid,ngrid)
      dimension px(ndec,2),pxm(ngrid,ngrid,ndec,2)
      common /grid/pxm,w
c
      do 5 i=1,ngrid
         w0(i) = 1.d0
         if((i.eq.1).or.(i.eq.ngrid)) w0(i) = 0.5d0
    5 continue
c
      dz=1.d0/float(ngrid-1)
      do 30 j=1,ngrid
         beta = dz*(j-1)
         do 20 k=1,ngrid
            if (k.lt.ngrid) then
               pz = 0.5d0 + 0.5d0*dz*(k-1)
            else
               pz = 0.999d0
            endif
            gam = 2.d0*dlog(pz/(1.d0-pz))
            call fpx(gam,beta,px)
            do 10 m=1,ndec
               pxm(j,k,m,1) = px(m,1)
               pxm(j,k,m,2) = px(m,2)
   10       continue
            w(j,k) = w0(j)*w0(k)
   20    continue
   30 continue
c
      return
      end
c***********************************************************************
      subroutine fpx(gam,beta,px)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12)
      dimension px(ndec,2)
c
      y1 = 10.d0
      y2 = 12.d0
      y3 = 15.d0
      y4 = 5.d0
      y5 = 6.d0
      y6 = 8.d0
      pjb = (1.d0+beta)/3.d0
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(1,1) = za/(za+zb)
      px(1,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(2,1) = za/(za+zb)
      px(2,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(3,1) = za/(za+zb)
      px(3,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(4,1) = za/(za+zb)
      px(4,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(5,1) = za/(za+zb)
      px(5,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(6,1) = za/(za+zb)
      px(6,2) = zb/(za+zb)
c
      za = y1*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(7,1) = za/(za+zb)
      px(7,2) = zb/(za+zb)
c
      za = y2*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(8,1) = za/(za+zb)
      px(8,2) = zb/(za+zb)
c
      za = y3*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(9,1) = za/(za+zb)
      px(9,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y4*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(10,1) = za/(za+zb)
      px(10,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y5*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(11,1) = za/(za+zb)
      px(11,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y6*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(12,1) = za/(za+zb)
      px(12,2) = zb/(za+zb)

c      do 10 k=1,ndec
c         write(*,'(i2,2g12.4)')k,(px(k,j),j=1,2)
c   10 continue
c
      return
      end
c***********************************************************************
      subroutine fch
      implicit real*8 (a-h,o-z)
      parameter (nplay=101,ndec=12,ngrid=51)
      dimension pxm(ngrid,ngrid,ndec,2),w(ngrid,ngrid)
      dimension pch(nplay,ngrid,ngrid)
      dimension is(nplay,ndec)
      common /grid/pxm,w
      common /mats/is
      common /pchm/pch

      do 50 id=1,nplay
         do 30 j=1,ngrid
            do 25 k=1,ngrid
               pin = 1.d0
               do 20 m=1,ndec
                  ix = is(id,m)
                  pin = pin*pxm(j,k,m,ix)
   20          continue
               pch(id,j,k) = pin
   25       continue
   30    continue
   50 continue
      return
      end
c***********************************************************************
      subroutine savephi
      implicit real*8 (a-h,o-z)
      parameter (np=7,ndec=12,ngrid=51,npts=ngrid*ngrid)
      dimension x(np),ph1(npts),ph2(npts),gl(npts)
      dimension pxm(ngrid,ngrid,ndec,2),w(ngrid,ngrid)
      common /grid/pxm,w
      common /glm/gl
      common /est/x
c
      b2  = 1.d0
      pz2 = 0.999d0
      sb2 = x(1)
      sz2 = x(2)
      b1  = x(3)
      pz1 = x(4)
      sb1 = x(5)
      sz1 = x(6)
      a2  = x(7)
      a1  = 1.d0 - a2
c
      gt2 = 0.d0
      dz = 1.d0/float(ngrid-1)
      do 20 j=1,ngrid
         j0 = ngrid*(j-1)
         bj = dz*(j-1)
         pb = (bj - b2)/sb2
         pb = pb**2
         pb = dexp(-0.5d0*pb)
         do 10 k=1,ngrid
            pk = 0.5d0 + 0.5d0*dz*(k-1)
            pg = (pk - pz2)/sz2
            pg = pg**2
            pg = dexp(-0.5d0*pg)
            ph2(j0+k) = w(j,k)*pg*pb
            gt2 = gt2 + ph2(j0+k)
   10    continue
   20 continue
c
      gt1 = 0.d0
      dz = 1.d0/float(ngrid-1)
      do 40 j=1,ngrid
         j0 = ngrid*(j-1)
         bj = dz*(j-1)
         pb = (bj - b1)/sb1
         pb = pb**2
         pb = dexp(-0.5d0*pb)
         do 30 k=1,ngrid
            pk = 0.5d0 + 0.5d0*dz*(k-1)
            pg = (pk - pz1)/sz1
            pg = pg**2
            pg = dexp(-0.5d0*pg)
            ph1(j0+k) = w(j,k)*pg*pb
            gt1 = gt1 + ph1(j0+k)
   30    continue
   40 continue
c
      gl(1) = a1*ph1(1)/gt1 + a2*ph2(1)/gt2
      do 50 i=2,npts
         gl(i) = gl(i-1) + a1*ph1(i)/gt1 + a2*ph2(i)/gt2
   50 continue
c      write(*,'(21f8.4)')(gl(i),i=1,npts)
c      stop
      return
      end
c***********************************************************************
      double precision function famoeb(x)
      implicit real*8 (a-h,o-z)
      parameter (np=7)
      dimension x(np)
c
      s0 = 0.05d0
      b2  = 1.d0
      pz2 = 0.999d0
      sb2 = 1.d0/(1.d0+dexp(x(1)))+s0
      sz2 = 9.95d0/(1.d0+dexp(x(2)))+s0
      b1  = 1.d0/(1.d0+dexp(x(3)))
      pz1 = 0.5d0/(1.d0+dexp(x(4)))+0.5d0
      sb1 = 1.d0/(1.d0+dexp(x(5)))+s0
      sz1 = 1.d0/(1.d0+dexp(x(6)))+s0
      a2  = 1.d0/(1.d0+dexp(x(7)))
c
      call fam(b2,pz2,sb2,sz2,b1,pz1,sb1,sz1,a2,pt)

      famoeb = -pt
c
      return
      end
c***********************************************************************
      subroutine fam(b2,pz2,sb2,sz2,b1,pz1,sb1,sz1,a2,pt)
      implicit real*8 (a-h,o-z)
      parameter (nplay=101,ndec=12,ngrid=51)
      dimension pch(nplay,ngrid,ngrid),phi1(ngrid,ngrid),
     &          phi2(ngrid,ngrid)
      common /pchm/pch

      call fprob(b2,pz2,sb2,sz2,b1,pz1,sb1,sz1,phi1,phi2)
      a1 = 1.d0-a2
c
      pt = 0.d0
      do 50 id=1,nplay
         pid = 0.d0
         do 20 j=1,ngrid
            do 15 k=1,ngrid
               phi = a1*phi1(j,k) + a2*phi2(j,k)
               pid = pid + phi*pch(id,j,k)
   15       continue
   20    continue
         if (pid .gt. 1.d-307) then
            pid = dlog(pid)
         else
            pid = -706.d0
         endif
         pt = pt + pid
   50 continue
      return
      end
c***********************************************************************
      subroutine fprob(b2,pz2,sb2,sz2,b1,pz1,sb1,sz1,phi1,phi2)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12,ngrid=51)
      dimension pxm(ngrid,ngrid,ndec,2),w(ngrid,ngrid)
      dimension phi1(ngrid,ngrid),phi2(ngrid,ngrid)
      common /grid/pxm,w
c
      dz = 1.d0/float(ngrid-1)
      pt2 = 0.d0
      do 20 j=1,ngrid
         bj = dz*(j-1)
         pb = (bj - b2)/sb2
         pb = pb**2
         pb = dexp(-0.5d0*pb)
         do 10 k=1,ngrid
            pk = 0.5d0 + 0.5d0*dz*(k-1)
            pg = (pk - pz2)/sz2
            pg = pg**2
            pg = dexp(-0.5d0*pg)
            phi2(j,k) = w(j,k)*pg*pb
            pt2 = pt2 + phi2(j,k)
   10    continue
   20 continue
      do 30 j=1,ngrid
         do 25 k=1,ngrid
            phi2(j,k) = phi2(j,k)/pt2
   25    continue
c      write(*,'(11f7.3)')(phi2(j,k),k=1,ngrid)
   30 continue
c
      pt1 = 0.d0
      do 50 j=1,ngrid
         bj = dz*(j-1)
         pb = (bj - b1)/sb1
         pb = pb**2
         pb = dexp(-0.5d0*pb)
         do 40 k=1,ngrid
            pk = 0.5d0 + 0.5d0*dz*(k-1)
            pg = (pk - pz1)/sz1
            pg = pg**2
            pg = dexp(-0.5d0*pg)
            phi1(j,k) = w(j,k)*pg*pb
            pt1 = pt1 + phi1(j,k)
   40    continue
   50 continue
      do 60 j=1,ngrid
         do 55 k=1,ngrid
            phi1(j,k) = phi1(j,k)/pt1
   55    continue
c      write(*,'(11f7.3)')(phi1(j,k),k=1,ngrid)
   60 continue
c
      return
      end
c***************************************************************
      subroutine data(it,dseed)
      implicit real*8 (a-h,o-z)
      parameter (nplay=101,ndec=12,ngrid=51,npts=ngrid*ngrid)
      dimension pxm(ngrid,ngrid,ndec,2),w(ngrid,ngrid)
      dimension gl(npts)
      dimension is(nplay,ndec)
      common /grid/pxm,w
      common /glm/gl
      common /mats/is
c
      do 50 id=1,nplay
         r = ggubfs(dseed)
         do 15 i=1,npts
            if (gl(i) .gt. r) goto 16
   15    continue
c
   16    j = (i-1)/ngrid + 1
         k = i - ngrid*(j-1)
c
         do 30 m=1,ndec
            x = ggubfs(dseed)
            p1 = pxm(j,k,m,1)
            if (x .le. p1) then 
               is(id,m)=1
            else
               is(id,m)=2
            endif
   30    continue
c      write(*,'(13i3)')id,(is(id,k),k=1,ndec)
   50 continue
      return
      end
c
c***********************************************************************
      subroutine init(xint,p,y,famoeb,delta)
      implicit real*8 (a-h,o-z) 
      parameter (np=7,mp=np+1)
      dimension q(np,np),p(mp,np),x(np),y(mp) 
      dimension xint(np),xx(np),xxx(np)
      external famoeb 
c 
      do 1 k=1,np
         xxx(k)=xint(k)
    1 continue
c
c      write(*,'(//''Initial Parameter Values:'')')
c      do 14 k=1,np
c         write(*,'(3x,g12.6,3x,g12.6)')xint(k),xxx(k)
c   14 continue
c
c      y0=-famoeb(xxx)
c      write(*,'(/''initial log-likelihood = '',g14.7/)')y0
c      stop
c
      p(1,1)=1.d0
      p(2,1)=-1.d0
      if (np.eq.1) goto 51
      do 5 k=2,np
         p(1,k)=0.d0
         p(2,k)=0.d0
    5 continue
c
      do 50 it=2,np
         f0=float(it)
         f1=(dsqrt(f0*f0-1))/f0
         f2=-1.d0/f0
c
         do 11 j=1,it
            do 10 k=1,np
               q(j,k)=p(j,k)
   10       continue
   11    continue
         do 15 k=1,np
            p(1,k)=0.d0
   15    continue
         p(1,it)=1.d0
         do 21 j=2,it+1
            p(j,it)=f2
            do 20 k=1,it-1
               p(j,k)=f1*q(j-1,k)
   20       continue
   21    continue
   50 continue
   51 continue
c
c Initial step size
c      delta=0.1d0
c      write(*,'(''  stepsize = '',g12.3//)')delta
c 
      do 80 j=1,mp
         do 79 k=1,np
            p(j,k) = xxx(k) + delta*p(j,k)
   79    continue
   80 continue
c 
c Initialize y. 
      do 92 i=1,mp 
         do 91 j=1,np 
            x(j)=p(i,j) 
   91       continue 
         y(i)=famoeb(x) 
c     write(*,'(i2,9(1x,g12.4),g14.7)')i,(x(j),j=1,np),y(i)
   92    continue 
      return 
      end 
c*********************************************************************
c ************** newlibrary of subroutines ***************************
      subroutine avg(p,x,np,mp)
      implicit real*8 (a-h,o-z)
      dimension p(mp,np),x(np)
c
c compute average values of vertices of the simplex.
      do 15 k=1,np
         do 14 j=2,mp
   14       p(1,k)=p(1,k)+p(j,k)
   15    x(k)=p(1,k)/mp
      return
      end
c------------------------------------------------------------------------
      subroutine amoeba(p,y,mp,np,ndim,ftol,funk,iter)
      implicit real*8 (a-h,o-z)
      parameter (nmax=20,alpha=1.d0,beta=0.5d0,gamma=2.d0)
      dimension p(mp,np),y(mp),pr(nmax),prr(nmax),pbar(nmax)
c
c This is the simplex algorithm from Numerical Recipes, with 
c modifications.
c When amoeba is called, iter must be set equal to the maximum
c number of iterations that will be allowed.  On return,
c iter equals the actual number of iterations.
c Thus, an example of a calling sequence would be:
c          itmax=5000
c          iter=5000
c          call aboeba(p,y,mp,np,ndim,ftol,funk,iter)
c          if (iter.ge.itmax) ...
c On return, if the if condition evaluates to true, convergence
c has not been acheived in the alloted number of iterations.
c     --PWW   1/26/94
c
c
      itmax=iter
      if (itmax.le.0) write(*,100)itmax
  100 format('*****WARNING: Maximum no. of iterations passed to',
     &       ' amoeba=',i6,'*****')
      mpts=ndim+1
      iter=0
    1 ilo=1
      if(y(1).gt.y(2))then
        ihi=1
        inhi=2
      else
        ihi=2
        inhi=1
      endif
      do 11 i=1,mpts
        if(y(i).lt.y(ilo)) ilo=i
        if(y(i).gt.y(ihi))then
        inhi=ihi
        ihi=i
        else if(y(i).gt.y(inhi))then
        if(i.ne.ihi) inhi=i
        endif
   11 continue
      rtol=2.d0*dabs(y(ihi)-y(ilo))/(dabs(y(ihi))+dabs(y(ilo)))
      if(rtol.lt.ftol)return
c
c check number of iterations.
      if (iter.ge.itmax) then
         write(*,101)itmax
         return
         end if
  101 format('*****Maximum iterations (',i6,') exceeded in amoeba.'/
     &       '     Control passed to calling routine.*****')
      iter=iter+1
      do 12 j=1,ndim
        pbar(j)=0.d0
   12 continue
      do 14 i=1,mpts
        if(i.ne.ihi)then
        do 13 j=1,ndim
        pbar(j)=pbar(j)+p(i,j)
   13   continue
        endif
   14 continue
      do 15 j=1,ndim
        pbar(j)=pbar(j)/ndim
        pr(j)=(1.d0+alpha)*pbar(j)-alpha*p(ihi,j)
   15 continue
      ypr=funk(pr)
      if(ypr.le.y(ilo))then
        do 16 j=1,ndim
        prr(j)=gamma*pr(j)+(1.d0-gamma)*pbar(j)
   16   continue
        yprr=funk(prr)
        if(yprr.lt.y(ilo))then
        do 17 j=1,ndim
        p(ihi,j)=prr(j)
   17   continue
        y(ihi)=yprr
        else
        do 18 j=1,ndim
        p(ihi,j)=pr(j)
   18   continue
        y(ihi)=ypr
        endif
      else if(ypr.ge.y(inhi))then
        if(ypr.lt.y(ihi))then
        do 19 j=1,ndim
        p(ihi,j)=pr(j)
   19   continue
        y(ihi)=ypr
        endif
        do 21 j=1,ndim
        prr(j)=beta*p(ihi,j)+(1.d0-beta)*pbar(j)
   21   continue
        yprr=funk(prr)
        if(yprr.lt.y(ihi))then
        do 22 j=1,ndim
        p(ihi,j)=prr(j)
   22   continue
        y(ihi)=yprr
        else
        do 24 i=1,mpts
        if(i.ne.ilo)then
        do 23 j=1,ndim
                pr(j)=0.5d0*(p(i,j)+p(ilo,j))
                p(i,j)=pr(j)
   23   continue
        y(i)=funk(pr)
        endif
   24   continue
        endif
      else
        do 25 j=1,ndim
        p(ihi,j)=pr(j)
   25   continue
        y(ihi)=ypr
      endif
      go to 1
      end
c**************************************************************
      subroutine sort(n,ra)
      implicit real*8 (a-h,o-z)
      dimension ra(n)
      l=n/2+1
      ir=n
10    continue
        if(l.gt.1)then
        l=l-1
        rra=ra(l)
        else
        rra=ra(ir)
        ra(ir)=ra(1)
        ir=ir-1
        if(ir.eq.1)then
        ra(1)=rra
        return
        endif
        endif
        i=l  
        j=l+l
20      if(j.le.ir)then
        if(j.lt.ir)then
        if(ra(j).lt.ra(j+1))j=j+1
        endif
        if(rra.lt.ra(j))then
        ra(i)=ra(j)
        i=j
        j=j+j
        else
        j=ir+1
        endif
        go to 20
        endif
        ra(i)=rra
      go to 10
      end      
c========================================================================
      double precision function ggubfs (dseed)
      double precision   dseed
      double precision   d2p31m,d2p31
      data              d2p31m/2147483647.d0/
      data              d2p31 /2147483648.d0/
      dseed = dmod(16807.d0*dseed,d2p31m)
      ggubfs = dseed / d2p31
      return
      end

